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Abstract 

Using an improved version of the previously introduced CRASH (Cosmic Ray Accel- 
eration SHock) code, we have calculated the time evolution of cosmic-ray (CR) mod- 
ified quasi-parallel plane shocks for Bohm-like diffusion, including self-consistent 
models of Alfven wave drift and dissipation, along with "thermal leakage injection" 
of CRs. The new simulations follow evolution of the CR distribution to much higher 
energies than our previous study, providing a better examination of evolutionary 
and asymptotic behaviors. The postshock CR pressure becomes constant after quick 
initial adjustment, since the evolution of the CR partial pressure expressed in terms 
of a momentum similarity variable is self-similar. The shock precursor, which scales 
as the diffusion length of the highest energy CRs, subsequently broadens approxi- 
mately linearly with time, independent of diffusion model, so long as CRs continue 
to be accelerated to ever-higher energies. This means the nonlinear shock structure 
can be described approximately in terms of the similarity variable, x/{ust), where Ug 
is the shock speed once the postshock pressure reaches an approximate time asymp- 
totic state. As before, the shock Mach number is the key parameter determining the 
evolution and the CR acceleration efficiency, although finite Alfven wave drift and 
wave energy dissipation in the shock precursor reduce the effective velocity change 
experienced by CRs, so reduce acceleration efficiency noticeably, thus, providing a 
second important parameter at low and moderate Mach numbers. For low Mach 
numbers (Mq ^ 5) the CR acceleration efficiency depends on the thermal leakage 
injection rate, the Alfvenic Mach number, and any preexisting CR population. How- 
ever, these dependences become weak for high shock Mach numbers of Mq > 30. 
To evaluate CR acceleration efficiencies in the simulated shocks we present for a 
wide range of shock parameters a "CR energy ratio", $(Mo), comparing the time 
asymptotic volume-integrated energy in CRs to the time-integrated kinetic energy 
flux through the shock. This ratio asymptotes to roughly 0.5 for sufficiently strong 
shocks. The postshock CR pressure is also approximately 1/2 the momentum flux 
through the shock for very high Mach numbers. 
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1 Introduction 

Astrophysical plasmas, from the interplanetary gas inside the heliosphere to 
the galaxy intracluster medium (ICM), are magnetized and turbulent and 
contain nonthermal particles in addition to gas thermal particles. So, under- 
standing complex interactions among these different components is critical to 
the study of many astrophysical problems. In coUisionless shocks entropy is 
generated via collective electromagnetic viscosities, i.e., interactions of charged 
particles with turbulent fields [33|. Some suprathermal particles of the shock 
heated gas can leak upstream, their streaming motions against the background 
fluid exciting MHD Alfven waves upstream of the shock |6|32j . Then those par- 
ticles can be further accelerated to very high energies through multiple shock 
crossings resulting from resonant scatterings with the self-excited Alfven waves 
in the flows converging across the shock [T3|10ll36j . 

Detailed nonlinear treatments of diffusive shock acceleration (DSA) account 
for incoming thermal particles injected into the CR population (e.g., P^37|21j ) 
as a consequence of incomplete thermalization by coUisionless dissipation pro- 
cesses. Those particles, while relatively few in number, can subsequently accu- 
mulate a major fraction of the shock kinetic energy as their individual energies 
increase p!6|26] . Such predictions are supported by a variety of observations 
including direct measurements of particle spectra at interplanetary shocks, 
nonthermal 7-ray, X-ray and radio emissions of supernova remnant shocks 
and also possibly the ICM of some X-ray clusters (e.g., p!0|l3|f¥T] ) . CR accel- 
eration may be universal to astrophysical shocks in diffuse, ionized media on 
all scales. 

Unlike an ideal gasdynamic shock, downstream states of a CR modified shock 
cannot be determined in a straightforward way by simple jump conditions 
across the shock front. This is because the shock transition depends on the 
CR pressure distribution upstream of the dissipative subshock. The particle ac- 
celeration takes place on diffusion time and length scales (td(p) = k{p)/u1 and 
ld{p) = ti{p)/us, where k{p) is diffusion coefficient and Ug is the shock speed), 
which are much larger than the shock dissipation scales. Unless or until some 
boundary condition limits the maximum CR momentum, this structure will 
continue to evolve along with the CR distribution. Thus the evolution of a CR 

Email addresses: kangOuju . es . pusan . ac . kr (Hyesung Kang), 
twj@astro.umn.edu (T. W. Jones). 

URL: www.astro.umn.edu/~twj (T. W. Jones). 



2 



shock with a finite age and size should properly be followed by time-dependent 
numerical simulations. In addition, complex interplay among CRs, resonant 
waves, and the underlying gas flow (i.e., thermal leakage injection, self-excited 
waves, resonant scatterings of particles by waves, and non-linear feedback to 
the gas flow) is model dependent and not yet understood completely. 

In the time dependent kinetic equation approach to numerical study of CR 
acceleration at shocks, the diffusion-convection equation for the particle mo- 
mentum distribution, f{p,x,t), is solved, along with suitably modifled gasdy- 
namic equations (e.g., [22] )• Since accurate solutions to this equation require 
a computational grid spacing smaller than the particle diffusion length, ld{p), 
and since realistic diffusion coefficients have steep momentum dependence, a 
wide range of length scales must be resolved in order to follow the CR accel- 
eration from the injection momentum (typically Pmj/'^pC ~ 10~^) to highly 
relativistic momenta {p/rripC ^ 1). This constitutes an extremely challenging 
numerical task, which can require rather extensive computational resources, 
especially if one allows temporal and spatial evolution of diffusion behavior. To 
overcome this numerical problem in a generally applicable way we have built 
the CRASH (Cosmic-Ray Acceleration SHock) code by implementing Adap- 
tive Mesh Refinement (AMR) techniques and subgrid shock tracking methods 
[25|30] in order to enhance computational efficiency. The CRASH code also 
treats thermal leakage injection self-consistently by adopting a shock trans- 
parency function for suprathermal particles in the shock |26j . 

We previously applied our CRASH code in a plane-parallel geometry to cal- 
culate the nonlinear evolution of CR modified shocks in the absence of sig- 
nificant local Alfven wave heating and advection [26|27ll29j . For those models 
the shock sonic Mach number, Mq, largely controlled the thermal leakage in- 
jection rate and the CR acceleration efficiency in evolving modified planar 
shocks, since Mq determines the relative velocity jump across the shock and 
consequently the degree of shock modification by CRs. In all but some of the 
highest Mach number shocks the CR injection rate and the postshock CR 
pressure approached time-asymptotic values when a balance was achieved be- 
tween acceleration/injection and diffusion/advection processes. This resulted 
in an approximate "self-similar" flow structure, in the sense that the shock 
structure broadened approximately linearly in time, so that the shock struc- 
ture could be expressed in terms of the similarity coordinate Ugt. It is likely 
that all of the models would have reached such asymptotic dynamical struc- 
tures eventually, but performance limits in the version of the code in use at 
that time prevented us from extending some of the simulations long enough to 
confirm that. The CR distribution evolved only to Pmax/'^pC ~ 10 — 10^. Based 
on the self-similar evolution reported in our previous work, we calculated the 
ratio of CR energy to infiowing kinetic energy, $, (see Eq. [12] below) as a 
measure of the CR acceleration efficiency in a time- asymptotic limit. The CR 
energy ratio, $, increased with the shock Mach number, but approached ~ 0.5 
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for large shock Mach numbers, Mq > 30, and it was relatively independent 
of other upstream properties or variation in the injection parameter. In those 
shocks where we observed time asymptotic dynamical behaviors the postshock 
CR pressures were ~30 - 60% of the ram pressure in the initial shock frame, 
this ratio increasing with Mach number. For some of the highest Mach num- 
ber shocks in that study CR pressure continued to increase to the end of the 
simulation, so the final values could not be measured. Finally, the presence of 
a preexisting, upstream CR population was seen in those earlier simulations 
to be equivalent to having slightly more efficient thermal leakage injection for 
such strong shocks, while it could substantially increase the overall CR energy 
in moderate strength shocks with Mq < 3. 

In the present paper, we revisit the problem of self-similar evolution of CR 
modified shocks with a substantially improved numerical scheme that enables 
us to follow the particle acceleration to energies much higher than we con- 
sidered before. This allows us to measure asymptotic dynamical properties 
for all the newly simulated shocks and to demonstrate that the self-similar 
evolution of the CR partial pressure in terms of a momentum similarity vari- 
able leads to the constancy of the postshock CR pressure. We also include in 
the new work the effects of Alfven wave drift and dissipation in the shocks. 
The time asymptotic CR acceleration efficiency is once again controlled by 
the shock Mach number but diminished as the ratio of Alfvenic Mach number 
to sonic Mach number decreases. The asymptotic shock properties are largely 
independent of the magnitude of the spatial diffusion coefficient and also its 
subrelativistic momentum dependence. 

The basic equations and details of the numerical method are described in §2. 
We present simulation results for a wide range of shock parameters in §3, 
followed by a summary in §4. 



2 Numerical Method 

2.1 Basic Equations 

The evolution of CR modified shocks depends on a coupling between the gas- 
dynamics and the diffusive CRs. That coupling takes place by way of resonant 
MHD waves, although it is customary to express the pondermotive wave force 
and dissipation in the plasma through the associated CR pressure distribu- 
tion properties along with a characteristic wave propagation speed (usually 
the Alfven speed) (e.g., [l2p] ). Consequently, in our simulations we solve the 
standard gasdynamic equations with CR pressure terms added in the conser- 
vative, Eulerian formulation for one dimensional plane-parallel geometry. The 
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evolution of a modified entropy, S = Pg/p^^~^, is followed everywhere except 
across the subshock, since for strongly shocked flows, numerical errors in com- 
puting the gas pressure from the total energy can lead to spurious entropy 
generation with standard methods, especially in the shock precursor |26j . 

dt + dx ^' 



where Pg and Pc are the gas and the CR pressure, respectively, Cg = Pg/[p(7g — 1)] + 
V? /2 is the total energy of the gas per unit mass. The remaining variables, ex- 
cept for L and W have standard meanings. The injection energy loss term, 
L{x,t), accounts for the energy carried by the suprathermal particles injected 
into the CR component at the subshock and is subtracted from the postshock 
gas immediately behind the subshock. Gas heating due to Alfven wave dissipa- 
tion in the upstream region is represented by the term W{x, t) = —VAdPc/dx, 
where va = B/y/Anp is the Alfven speed. This commonly used dissipation 
expression derives from a quasi-linear model in which Alfven waves are ampli- 
fied by streaming CRs and dissipated locally as heat in the precursor region 
(e.g., [22]). 

The CR population is evolved by solving the diffusion-convection equation in 
the form, 

, .dg 19, ..dg ^ , d . , .dg. 



where g = p'^f, with f{p,x,t) the pitch angle averaged CR distribution, and 
where y = ln(p), while k{x,p) is the spatial diffusion coefficient [12]. For sim- 
plicity we always express the particle momentum, p in units rripC and consider 
only the proton CR component. The wave speed is set to be = va in the 
upstream region, while we use = in the downstream region. This term 
reflects the fact that the scattering by Alfven waves tends to isotropize the CR 
distribution in the wave frame rather than the bulk- flow, gas frame [12]. Up- 
stream, the waves are expected to be dominated by the streaming instability, 
so face upwind. Behind the shock, various processes, including wave reflection, 
are expected to lead to a more nearly isotropic wave field (e.g., [2]). 
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Eqs. (l)-(5) are simultaneously integrated by the CRASH code in plane- 
parallel geometry. The detailed numerical description can be found in Kang 
et al. 2002 ||26j. A key performance feature of the CRASH code is multiple 
levels of refined grids (typically Ig = 8 — 10) strategically laid around the sub- 
shock to resolve the diffusion length scale of the lowest energy particles near 
injection momenta. Grid refinement spans a region around the subshock just 
large enough to include comfortably the diffusion scales of dynamically impor- 
tant high energy CRs with enough levels to follow freshly injected low energy 
CRs with sufficient resolution to produce converged evolutionary behaviors. 
To accomplish grid refinement effectively it is necessary to locate the subshock 
position exactly. Thus, we track the subshock as a moving, discontinuous jump 
inside the initial, uniform and fixed grid [25] . 



2.2 Diffusion Model 



We considered in this study two common choices for diffusion models. First 
is the Bohm diffusion model, which represents scattering expected for a satu- 
rated wave spectrum and gives what is generally assumed to be the minimum 
diffusion coefficient as kb = l/3rgi;, when the particles scatter within a path 
of one gyration radius (i.e., Amfp ~ Tg). This gives 

^^^P^ = ^^^ (^2 + 1)1/2 - (6) 



The coefficient Kn = mc^/{3eB) = (3.13 x 10^^011^8-^)3-^, where 5^ is the 
magnetic field strength in units of microgauss. There has been much discussion 
recently about amplification of the large scale magnetic field within the shock 
precursor (e.g., [52]|35|44] ). Since physical models of that evolution are still 
not well developed, we will assume for simplicity in the simulations presented 
here that the large scale field is constant through the shock structure. 

Because of its steep momentum dependence in the nonrelativistic regime, the 
Bohm diffusion model requires an extremely fine spatial grid resolution when- 
ever nonrelativistic CRs are present. On the other hand the form of K,{p) for 
nonrelativistic momenta mostly impacts only early evolution of CR-modified 
shocks, when CR feedback is dominated by nonrelativistic particles and ther- 
mal leakage injection rates are adjusting rapidly to changes from initial condi- 
tions. So, to concentrate computational effort more efficiently, we adopted in 
some previous works [25)130] a "Bohm-like" diffusion coefficient that includes 
a weaker momentum dependence for the non-relativistic regime, 

«bl(p) = «nP- (7) 
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According to those previous studies, the differences in results between the two 
models are minor except during early nonlinear shock evolution, as expected. 
Thanks to the weaker momentum dependence of kbl we can, for given com- 
putational resources, calculate numerically converged models with smaller 
resulting in the acceleration of CRs to higher momenta. 

In order to quench the well-known CR acoustic instability in the precursor of 
highly modified CR shocks (e.g., [21]), we assume a density dependence for 
the diffusion coefficient, {p/po)~^, so that k{x,p) = kb{p/po)~^ or n{x,p) = 
kbl(p/Po) S where po is the upstream gas density. This density dependence 
also models enhancement of the Alfven wave magnetic field amplitude due to 
flow compression. We note, also, for clarity that hereafter we use the subscripts 
'0', '1', and '2' to denote conditions far upstream of the shock, immediately 
upstream of the gas subshock and immediately downstream of the subshock, 
respectively. 



2.3 Thermal Leakage Model 

In the CRASH code suprathermal particles are injected as CRs self-consistently 
via "thermal leakage" through the lowest CR momentum boundary. The ther- 
mal leakage injection model emulates the filtering process by which suprather- 
mal particles well into the tail of the postshock Maxwellian distribution leak 
upstream across the subshock [371136] . This filtering is managed numerically 
by adopting a "transparency function", resc(e_B, t;), that expresses the proba- 
bility of supra-thermal particles at a given velocity, v, successfully swimming 
upstream across the subshock through the postshock MHD waves |2TP5] . The 
one model parameter, = Bo/B±, is the ratio of the amplitude of the post- 
shock wave field interacting with the low energy particles, B±, to the general 
magnetic field, Bq, which is aligned with the shock normal in these simulations. 
The transparency function fixes the lowest momentum of the CR component 
in our simulations from the condition that Tesc > (i.e., non-zero probability 
to cross the subshock for CRs) for p > pi, where pi = (u2/c)(l + eB)/€B and 
U2 is the downstream flow speed in the subshock rest frame. Initially pi is 
determined by the downstream speed of the initial shock, but it decreases as 
the subshock weakens and then it becomes constant after the CR modified 
shock structures reach asymptotic states. 

Since suprathermal particles have to swim against the scattering waves ad- 
vecting downstream, the subshock Mach number, Mg, is one of the key shock 
characteristics that control the injection fraction. Previous simulations showed 
that injection is less efficient for weaker shocks, but becomes independent of 
Mq for strong shocks, when the subshock compression asymptotes [26]. For a 
given total shock Mach number, Mq, on the other hand, the injection rate 
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is controlled mainly by the parameter In practice we have found that 
~ 0.2 — 0.25 leads to an injection fraction in the range ~ 10~^ — 10~^. 
This is similar to commonly adopted values in other models that employ a 
fixed injection rate (e.g., pi5^IH] ). Although somewhat higher field turbulence 
values (0.25 < < 0.3) are suggested theoretically for strong shocks [M] . 
these evoke a start-up problem in the numerical simulations, since they lead 
to very rapid initial injection that cools the postshock flow too strongly for 
it to remain numerically stable. Once the shock structure becomes nonlinear, 
however, those influences moderate greatly, so that, as we have shown pre- 
viously, the ultimate CR acceleration behavior depends only weakly on e^. 
In fact, we have found previously for strong shocks that the time-asymptotic 
behaviors are very weakly dependent on [25], so its chosen value will have 
no influence on our final conclusions. 

We directly track the fraction of particles injected into the CR population as 
follows: 



where /cr is the CR distribution function above pi, while noUsflt is the number 
of particles passed through the shock until the time t. The highest momen- 
tum of the CR component, p2, is chosen so that it is well above pmax at the 
simulation termination time, where Pmax is defined in §2.4. 

2.4 CR Acceleration Efficiency 

The postshock thermal energy in a gasdynamic shock can, of course, be cal- 
culated analytically by the Rankine-Hugoniot jump condition. On the other 
hand, the CR population and the associated acceleration efficiency at CR 
modified shocks should properly be obtained through time-dependent integra- 
tion of the shock structure from given initial states, since the CR distribution 
depends on the shock structure, which is not discontinuous and continues to 
evolve so long as the CR population evolves. In particular CR modified shocks 
contain a smooth precursor in the upstream region whose scale height grows 
in time in proportion to the diffusion length of energetically dominant par- 
ticles. The total shock compression may, similarly, evolve over a significant 
time period. The standard expression for the mean acceleration timescale for 
a particle to reach momentum p in the test-particle limit of DSA theory is 
given by [31] 



/dxffM7r/cR(p,x,tydp 

UoUsfit 




Taccip) 
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In the test particle limit the shock compression is given by the Rankin- 
Hugoniot condition. Then assuming a 7g = 5/3 gasdynamic shock and a 
diffusion coefficient taking the density dependence indicated at the end of 
§2.2, this leads to 

r (p)^8 ^° ^ (10) 
racc[P) ~ » 2 _ 1 ^2 ' l^uj 



where Ug and Mq are the shock speed and sonic Mach number, respectively. 
While this expression should strictly speaking be modified in highly modified 
CR shocks, since the shock structure and the associated CR transport are more 
complex than assumed in Eq. ffTOj) ^llj , we find it empirically to be reasonably 
consistent with our results described below. Accordingly, we may expect and 
confirm below that the time- dependent evolution of our CR modified shocks 
will be determined primarily by the shock Mach number and can be expressed 
simply in terms of diffusion length and time scales. 

Within this model the highest momentum expected to be accelerated in strong 
shocks by the time t is set according to Eq. flTUl) by the relation t ~ 8K(pmax)/w^. 
In that case the scale height of the precursor or shock transition structure 
grows linearly with time as 

, 't(Pmax) 1 , 

tshock ~ ~ Q^st, I J- J- J 



independent of the magnitude or the momentum dependence of This 
evolution should continue until some other physics limits the increase in CR 
momentum, such as the finite size of the shock system. Since the CR pressure 
approaches a time-asymptotic value (see Figs. [SJIS] below), the evolution of 
the CR-modified shock becomes, under these circumstances, approximately 
self-similar, independent of the form of the diffusion coefficient, while 4hock 
grows linearly with time [27|28|29j . On the other hand, ^(pmax) ~ 4hockWs, 
so for Bohm-like diffusion, Pmax ~ 4hock(^)^s/^n- Thus, at a given time the 
CR distribution, g{p), extends for Bohm-like diffusion according to Praax ~ 

ult/{SKn) OC 

Fig. [1] shows a comparison of three models of a Mq = 20 shock with ks = 
KnP^ / + 1 using kn = 0.1, and with kbl = knP using k„ = 10"^ and 10~^ 
in units defined in the following section. The upper left panel, displaying the 
evolution of postshock CR pressure, demonstrates similar time-asymptotic 
values for all three models. At early times the CR pressure evolution depends 
on details of the model, including numerical properties such as spatial and 
momentum grid resolutions, and the previously described injection suppression 
scheme used to prevent start up problems. 
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The other panels in Fig. [T] illustrate shock structure comparisons for the three 
models at the end of the simulations. The shock structures, p{x) and Pc{x), 
are very similar, while the CR spectrum extends to different values of Pmax, 
inversely proportional to 

The self-similar evolution of CR modified shocks makes it useful to apply the 
ratio of CR energy to a fiducial kinetic energy flux through the shock as a 
simple measure of acceleration efficiency; namely. 



More specifically, this compares the total CR energy within the simulation box 
to the kinetic energy in the initial shock frame that has crossed the shock at a 
given time. As the shock structures approach time-asymptotic forms the above 
discussion suggests that ^(t) also may approach time-asymptotic values. This 
is confirmed in our simulations. We see also that the asymptotic $ values 
depend in our simulations primarily on shock sonic Mach number and are 
independent of n. 

The highest momenta achieved in our simulations are set by practical limits on 
computation time controlled by the vast range of diffusion times and lengths 
to be modeled. Still, the asymptotic acceleration efficiency ratio is almost 
independent of the maximum momentum reached. For the three models shown 
in Fig. [H for example, Pmax ~ 10, 10^, and 10^ at t = 10, depending on the 
value of Kn, but the CR energy ratio approaches similar values of $ ~ 0.4 for 
all three models. 



3 Simulation Set Up and Model Parameters 

3.1 Units and Initial Conditions 

The expected evolutionary behaviors described above naturally suggest conve- 
nient units. For example, given a suitable velocity scale, u, and length scale, x, 
a time scale t = x/u and diffusion coefficient scale, k = xu = u^t are implied. 
Alternatively, one can select the velocity scale along with a convenient scale 
for the diffusion coefficient, k, leading to a natural length scale, x = k/u and 
a related time scale, t = k/u^ = x/u. We will follow the latter convention in 
our discussion. In addition, given an arbitrary mass unit, m, which we take 
to be the proton mass, we can similarly normalize mass density in terms of 
f> = m/x^. Pressure is then expressed relative to pu^. For clarity we henceforth 




0.5poMs,0' 



(12) 
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indicate quantities normalized by the above scales using a tilde; for example, 
u, t, and kn- 

We start each simulation with a pure gasdynamic, right-facing shock at rest 
in the computational grid. We use the upstream gas speed in this frame, uq 
as our velocity scale, so that the initial, normalized shock speed is Usfi = 1 
with respect to the upstream gas. The upstream gas is specified with one 
of two temperature values, either Tq = lO^K or To = 10^, which represent 
warm or hot phase of astrophysical diffuse media, respectively. In astrophysical 
environments, for example, photoionized gas of 10^ K is quite common. Hot 
and ionized gas of > 10^ K is also found in the hot phase of the ISM [38P3] . The 
shock speed and the upstream temperature are related through the sonic Mach 
number, Mq, by the usual relation Usfl = Cs^Mq = 15 km s~"'^(To/10'^)"'^/^Mo = 
Mo, where Cg^ is the sound speed of the upstream gas. So, by choosing To and 
Mo, we set the physical value of the shock speed, which, in turn, determines the 
postshock thermal behavior. For CR distribution properties it is also necessary 
to define the speed of light, c, in terms of Pk = uo/c. The normalized upstream 
gas density is set to unity; i.e., po = 1. 

The preshock pressure is determined by the shock Mach number, Pgo = 
(l/7g)Mo~^, where the gas adiabatic index, 7g = 5/3. The postshock states 
for the initial shocks are determined by the Rankine-Hugoniot shock jump 
condition. For models with To = IC^K, Mo =10-80 is considered, since the gas 
would not be fully ionized at slower shocks {us < 150 km s~^), the postshock 
gas would often become radiative and the CR acceleration become inefficient 
owing to wave dissipation from ion-neutral collisions (e.g., [15j). For models 
with To = lO^K Mo = 2 - 30 (300 km s"^ < < 4500 km s'^) is considered, 
since the CR acceleration should be relatively independent of To for shocks 
with Mo > 30. 

In order to explore effects of pre-existing CRs, we also consider, as we did 
in our earlier work, models with To = 10® K that include an ambient (up- 
stream) CR population, f{p) oc for pi < p < P2 and set its pressure 
Pc,o = (0.25 — 0.3)Pgo. For these models, we adopt kb = 0.1p^/a/p^~-FT, 
Pi = (^s,o/c)(l + ^b)/^b and p2 = 10^. For strong shocks the presence of pre- 
existing CRs is similar in effect to having a slightly higher injection rate, so the 
time asymptotic shock structure and CR acceleration efficiency depend only 
weakly on such a pre-existing CR population [26||28] . On the other hand, for 
weak shocks, a pre-existing CR pressure comparable to the upstream gas pres- 
sure represents a significant fraction of the total energy entering the shock, 
so pre-existing CRs obviously have far more impact. In addition, the time 
asymptotic CR acceleration efficiency in weak shocks depends sensitively on 
the injection rate, so increases with increased shock transparency, controlled 
through cb (see §2.3). Hence, as we found in [29], we expect relatively weak 
CR shocks (Mo < 5) to be substantially altered by the presence of a finite 
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upstream Pc. 



3.2 Wave Drift and Heating 



As shown in earlier works (e.g., [2211^ and references cited therein), the CR 
acceleration becomes less efficient when Alfven wave drift and heating terms 
are included in the simulations. This behavior comes from two effects pre- 
viously mentioned in §2.1, both of which derive from the resonance inter- 
action between CRs and Alfven waves in the shock precursor. The Alfven 
waves stimulated by CR streaming in the precursor will propagate in the up- 
stream direction, so that the effective advection speed of the CRs into the 
subshock is reduced. In addition, if the energy extracted from CRs to am- 
plify these waves is locally dissipated, the heating rate in the precursor is 
increased with respect to the adiabatic rate, so that gas entering the sub- 
shock is relatively hotter, and the subshock strength is accordingly reduced. 
The significance of the effects depends on the sonic Mach number, Mq, rel- 
ative to the shock Alfvenic Mach number. Ma = Us/va] i-e., on the ratio of 
the Alfven speed to the sound speed (see §3.4). In a parallel shock we can 

write Mq/Ma = va/cs = \J26/ [7g(7g — 1)], where we introduce a convenient 
"Alfven parameter" as follows, 

» = = (13) 



This expresses the relative upstream Alfven and sound speeds in terms of 
the magnetic to thermal energy density ratio. The parameter f]p is the usual 
"plasma /?" parameter. For 7g = 5/3, /3p = 2/(39). and 9 = (5/9)(Mo/Ma)2. 
We emphasize for clarity that the present simulations are of parallel shocks, 
so that the direct dynamical role of the magnetic field has been neglected. 
Observed or estimated values are typically 6* ~ 0.1 for intracluster media and 
~ 1 for the interstellar medium of our Galaxy (e.g., [5IT2] ). So we consider 
0.1 < ^ < 1, and we will provide comparison to the weak field limit, ^ = 0. 
Fig. [2] concisely illustrates the importance of Alfven wave drift and heating 
effects on a Mach 10 shock. One can see that the postshock CR pressure, 
for example, decreases by more than a factor two when the sonic and Alfven 
Mach numbers become comparable. Most of the model results we show here 
used 9 = 0.1. For the Mach 10 shocks illustrated in Fig. [2]the associated wave 
terms have reduced the asymptotic CR pressure by about 30% with 6* = 0.1 
compared to the shock with no such terms included. 
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3.3 Grid Resolution and Convergence 



According to our previous studies [23j, the spatial grid resolution should be 
much finer than the diffusion length of lowest energy particles near the injec- 
tion momenta, i.e., Ax < 0.1ld{pinj)- Kang & Jones (2006), however, showed 
that the spherical, comoving CRASH code achieve good numerical conver- 
gence, even when Ax > ld{Pmi)- This was due to the fact in the execution 
of that code the gas subshock remains consistently inside the same comoving 
grid zone. In order to gain this benefit for our present simulations, we have 
modified our plane-parallel CRASH code so that the shock is again forced to 
remain inside the same refined grid zone by regularly redefining the under- 
lying Eulerian grid. The simulations employ eight levels of refinement with 
the grid spacing reduced by an integer factor of two between refinement lev- 
els. The spatial grid resolution on the coarsest, base grid is Axq = 2 x 10"'^, 
while on the finest, S*'', grid Axg = 7.8 x 10"^ With k„ = 10"^ the diffu- 
sion length for injection momenta pi ~ 10~^ becomes Id ~ 10^* for models 
with Bohm-like diffusion, kbl- Although Axg > ld{Pmi), we confirmed that the 
new plane-parallel CRASH code also achieves good numerical convergence in 
the simulations presented here. This improvement enables us to extend these 
simulations to CR momenta several orders of magnitude greater that those 
discussed in [29] . 

When solving the diffusion-convection equation, we used 230 - 280 uniformly 
spaced logarithmic momentum bins in the interval y = Inp = [Inpi, lnp2]- 



3.4 Results 



We show in Fig. [3] the time evolution of a Mq = 10 shock with Tq = lO^K, 
kbl = 10~^p, and 6 = 0.1. The wave amplitude parameter in the thermal leak- 
age model was assumed to be = 0.2, unless stated otherwise. The lower left 
panel follows evolution of the volume integrated CR distribution function rel- 
ative to the total number of particles (mostly in the thermal population) that 
have passed through the shock, i.e., G{p)/{nQUsfit), where G{p) = J dxg{x,p). 
As the CR pressure increases in the precursor in response to thermal leakage in- 
jection at the subshock and subsequent Fermi acceleration, the subshock weak- 
ens. The injection process is self-regulated in such a way that the injection rate 
reaches and stays at a nearly stable value after quick initial adjustment. Con- 
sequently, the postshock CR pressure reaches an approximate time-asymptotic 
value once a balance is established between fresh injection/acceleration and 
advection/diffusion of the CR particles away from the shock. 
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The CR pressure is calculated as 

Pc = Y"^pC^ J g{p)-^^=dlnp, (14) 
Pi 

so we define D{p) = g{p)p/ \/p^ + 1 as a 'partial pressure function'. The upper 
left panel of Fig. H] shows the evolution of D{p,Xs) at the subshock for the 
model shown in Fig. [31 Since D{p) stretches self-similarly in momentum space, 
we define a new momentum similarity variable as 

V - Hp/Pi) ^^5^ 

ln[Pmax(t)M]' 

The lowest momentum pi becomes constant after the subshock structure be- 
comes steady. The momentum at which numerical values of D{p) peaks is 
chosen as Pma^it) and similar to what is estimated approximately by applying 
the test-particle theory in §2.4. For the model shown in Fig. 3, pi ~ 0.01 and 
p^ax ~ 7.27 X 10^ i. 

We then define another 'partial pressure function', 

F{Z) = g{Z)-j^==ln[p^Ut)/p^] = DiZ)\n[p^Ut)/Pi]- (16) 

Its time evolution is shown as a function of Z at the upper right panel of Fig. 
m The plot demonstrates that the evolution of F{Z) becomes self-similar for 
i > 2. Since Pc oc /^^ D{p)d\np oc jQF{Z)dZ, the areas under the curves of 
D{p) or F{Z) in Fig. 4 represent the CR pressure at the shock. So the self- 
similar evolution of F{Z) implies the constancy of Pc,2- In this case Pc,2 ~ 0.31 
for t > 2. /^From that time forward the spatial distribution of Pc expands 
approximately linearly with time, as anticipated in §2.4. This demonstrates 
that the growth of a precursor and the shock structure proceed approximately 
in a self-similar way once the postshock CR pressure becomes constant. In the 
lower panels of Fig. 4 we also show the evolution of D{p) and F{Z) for another 
model with Mq = 50, Tq = lO^K, kb = OMp/y/jFTT and ^ = 0.1. For the 
this model pi ~ 2.0 x 10^'^ and Pmax ~ lOt. In this case, the CR pressure is 
mostly dominated by relativistic particles. 

Fig. O compares the CR distributions at t = 10 for two sets of models spanning 
a range of sonic Mach numbers for different diffusion model choices. For the 
simulations represented on the left, the diffusion coefficient is Bohm-like, with 
kbl = 10~^j9 and Tq = lO^K. The results on the right come from a Bohm 
diffusion model with ks = O.Olp^/ \/p^ + 1 and Tq = 10^ K. For all models in 
this figure, d = 0.1 and e = 0.2. The top panels show the CR distributions 
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at the shock, g{p,Xs) = p^f{p,Xs), while the middle panels show the spatially 
integrated G{p). The slopes of the integrated spectra, q = —d(\nGp)/d\np + 
4, are shown in the bottom panels. For strong shocks with Mq > 10, both 
p'^fip, Xg) and G{p) exhibit concave curvature at high momentum familiar from 
previous studies e.g., [9ll36lH] . The hardened high momentum slopes reflect the 
fact that higher momentum CRs have longer mean free scattering paths, so 
encounter an increased compression across the shock precursor and a greater 
velocity jump. 

It is interesting to note in Fig. [S] that the spatially integrated spectra in the 
stronger shocks show more obvious hardening at high momenta than do the 
spectra measured at the subshock. This is another consequence of the fact that 
ld{p) increases with momentum, so that the CR spectrum hardens considerably 
as one measures it further upstream of the subshock. CRs escaping upstream 
from such a shock would be dominated by the highest momentum particles 
available. 

Fig. [6] compares the evolution of shock properties for the same set of mod- 
els whose CR spectra are shown in Fig. O The evolution of density increase 
through the precursor, cXp = Pi/po, and the total compression, cr^ = P2/P0) 
are shown in each top panel. Compression through the subshock itself can be 
found through the ratio a.^ = P2/ Pi- We will discuss shock compression results 
in more detail below. 

The middle panels of Fig. O show evolution of the postshock gas and CR pres- 
sures normalized to the ram pressure at the initial shock. After fast evolution 
at the start, these ratios approach constant values for t > 1. The bottom 
panels monitor the thermal leakage injection fraction, ^{t). With the adopted 
value oieB = 0.2, the time asymptotic value of this fraction is ^ ~ lO"'* — 10~^ 
with the higher values for stronger shocks. We note that the behavior of ^ 
during the early phase (t < 1) is controlled mainly by the injection suppres- 
sion scheme used to prevent start up problems. If the diffusion of particles 
with lowest momenta near pi is better resolved with a finer grid spacing, one 
may observe some initial reduction of ^{t) as the subshock weakens in time 
(e.g.. Fig. 6 of The various plots show that postshock properties approach 
time-asymptotic states that depend on the shock Mach number. Generally, as 
is well known, the shock compression, the normalized postshock CR pressure 
and the thermal leakage injection fraction increase with shock Mach number. 
Nonlinear effects, illustrated here by increased shock compression and dimin- 
ished postshock gas pressure, are relatively unimportant for Mach numbers 
less than 10 or so. At high Mach number upstream gas temperature is rela- 
tively unimportant. In addition, noting that the simulations in the left and 
right panels employed different models for the low momentum diffusion co- 
efficient, the similarity between the analogous left and right plots illustrate 
the point made earlier that asymptotic dynamical behaviors are insensitive to 
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details in the diffusion coefficient. 



The Mach number dependences of several important shock dynamical prop- 
erties are illustrated in Fig. [71 These include the time-asymptotic values of 
at, (Tp, and the postshock CR and gas pressures, both normalized by the ini- 
tial momentum ffux, pou^Q. Dotted lines provide for comparison the 7g = 5/3 
gasdynamic compression ratio and postshock pressure. The strong shock lim- 
its of the gasdynamic compression ratio and normalized postshock pressure 
are 4 and 0.75, respectively. The two sets of models shown in Figs. 5-6 are 
shown, along with an additional set of models with To = lO^K, kbl = lO^^p, 
6 = 0.1, and = 0.25. The normalized postshock CR pressure increases 
with Mq, asymptoting towards Pc,2 ~ 0.5 at high Mach numbers. These val- 
ues are reduced ~ 30% from the analogous results presented in [29j, which 
is consistent with the comparison for the Mach 10 shock shown in Fig. [2l 
The normalized gas pressure simultaneously drops in response to increases 
in Pc, as one would anticipate. It is notable that it actually falls well below 
0.25, where one might naively expect it to asymptote in order to maintain 
a constant total postshock pressure with respect to pou^o- Instead it falls as 
Pg,2 = Pg,2/{poulo) ~ 0.4(Mo/10)-°■^ without flattening at high Mach num- 
bers. As we shall see below, this trend is consistent with expected evolution 
of gas in the precursor. Thus the total postshock pressure is less than that of 
a gasdynamic shock of the same initial Mach number, Mq, despite a softening 
of the equation of state coming from the CRs. 

For the = 0.2 cases shown in Fig. [7] the precursor, subshock, and to- 
tal compression ratios can be approximated by ~ 1.5(Mo/10)°'^^, as ~ 
3.3(Mo/10)°•°^ and at = apas ~ 4.9(Mo/10)°-33, respectively, for Mq > 5. 
For these models we also find the subshock Mach number ranges 3 < Mi(= 
Wi/cs^i) < 4 and depends on total Mach number roughly as Mi oc Mq-^. Since 
the Rankine-Hugoniot-derived gas shock compression in this Mach number 
range scales as ag oc M°'^^, so as oc (Mq'^)"'^^ is consistent with the above re- 
lation. Although weak, the variation of the subshock strength with full shock 
sonic Mach number is important to an understanding of the simulated flow be- 
haviors seen in our shock precursors. We note that the subshock behavior seen 
here in highly modified strong shocks, particularly that it generally evolves to- 
wards Mach numbers close to three, is consistent with previous analytic and 
numerical results, although in some other studies the subshock strength is 
completely independent of Mq (e.g., [7|l9|l36] ). The subshock strength in our 
simulations is determined by complex nonlinear feedback involving the ther- 
mal leakage injection process. As noted for the results in Fig. the injection 
rate, ^, generally increases with Mq. This simultaneously increases Pc, pro- 
viding extra Alfvenic heating in the precursor, while cooling gas entering the 
subshock, as energy is transferred to low energy CRs. It should not be surpris- 
ing that the balance of these feedback processes is not entirely independent 
of total Mach number. We comment in passing that the compression values 
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shown in Fig. [7] depend slightly on in the sense that larger values of 
result in a little higher CR injection rate and greater P^- 



Compression through the shock precursor has been discussed by several previ- 
ous authors (e.g., [7|H|33f26j ). We outline the essential physics to facilitate an 
understanding of our results. In a steady flow {d/dt = 0) the modified entropy 
equation (jl]) can be integrated across the precursor with 7^ = 5/3 to give 
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The same relation can be derived from a Lagrangian perspective applying the 
second law of thermodynamics to a parcel of gas flowing through the precur- 
sor, not necessarily in a steady state. For these simulations \W\ = \vAdPc/dx\. 
Given the previously mentioned result, Pc,2 = Pc,i ~ O.Spo'^s we can estimate 
that I ~ va/us ~ 1.34^/6 /Mq. Eq. f|T7|) is then similar to an expression given 
in |7]. When Alfven wave dissipation is small, so that 9 « 5/(4Mq), Eq. f|T71) 
gives the result dp ~ Mq^'^, appropriate for adiabatic compression and con- 
sistent with behaviors found in a number of previous analytic and numerical 
studies (e.g., [71l8]l33ll26j ). For the simulations represented in Figs.[3H7]the oppo- 
site limit actually applies, since ^ = 0.1 > 5/(4Mq) whenever Mq > 3.5. Then 
the strong Alfven dissipation limit of Eq. ([TT]) predicts cxp ~ 
Substituting our observed relation between subshock Mach number and to- 
tal Mach number. Mi oc Mq'^, we establish for fixed 6 an expected behavior, 
cTp oc Mq'^, very close to what is observed. If instead of 6 we had parame- 
terized the Alfven wave influence in terms of the Alfvenic Mach number, the 
analogous behavior of Eq. (1171) would have been ap oc 1/{M^J^M^^^). Then for 
fixed Alfvenic Mach number the precursor compression would vary only with 
the subshock Mach number, which, has only a weak dependence on the full 
shock Mach number. That agrees with results presented by [7j, for example. 

We mention that the precursor/subshock compression properties also explain 
the observed inverse dependence of postshock gas pressure on initial Mach 
number. In fact, in a steady flow it is easy to show that 
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Thus we expect an inverse relation between precursor compression and nor- 



17 



malized postshock gas pressure, close to what is observed in our simulations. 

Finally, Fig. [8] shows time-asymptotic values of the CR energy ratio, $(Mo), 
for models with different 9 (top panel) and different and pre-existing CRs 
(bottom panel). Models are shown with both Tg = lO^K and Tq = lO'^K. This 
figure demonstrates that the CR acceleration depends on the specific model 
parameters considered here. For models with ^ = 0.1 and T = 10^ K, the 
acceleration efficiency is reduced by up to ~ 50 % in comparison to models 
without Alfven wave drift and dissipation {6 = 0). For larger Alfven speeds 
with 6* ~ 1, the acceleration efficiency is reduced even more significantly as a 
consequence of strong preshock Alfvenic heating and wave advection. On the 
other hand, for models with Tq = lO^K and Mq > 20 the reduction factor is 
less than 15 % for 6* ~ 0.5. As shown in previous studies [26|27] . larger values 
of es lead to higher thermal leakage injection and so more CR energy. Also 
the presence of preexisting CRs facilitates thermal leakage injection, leading 
to more injected CR particles and higher acceleration efficiency. Thus, an 
accurate estimate of the CR energy generated at quasi-parallel shock requires 
detail knowledge of complex physical processes involved. Fortunately, all theses 
dependences become gradually weaker at higher Mach numbers, and $ tends 
to approach 0.5 for Mq > 30. It seems likely that this asymptotic efficiency 
would also apply for sufficiently large 6. For low Mach number shocks, it is 
not yet possible to make simple, model-independent efficiency predictions. 



4 Conclusion and Summary 

While full thermalization takes place instantaneously at a simple, discontinu- 
ous jump in an ideal gasdynamic shock, CR acceleration and the corresponding 
modifications to the underlying flow depend on suprathermal particles passing 
back and forth diffusively across the shock structure via resonant scatterings 
with MHD waves. So the time dependent evolution of CR modified shocks 
depend on complex interactions between the particles, waves in the magnetic 
field, and underlying plasma flow. These processes develop on the diffusion 
time scale and diffusion length scale, which are expected to be increasing 
functions of particle momentum, so CR acceleration and shock evolution rates 
slow over time. Here we have carried out time-dependent DSA simulations to 
study the evolution of plane CR modified shocks with quasi-parallel magnetic 
fields. Simple models for thermal leakage injection and Alfven wave propaga- 
tion and dissipation are included. By adopting an underlying base grid that 
moves with the subshock, our CRASH code with shock tracking and AMR 
techniques can achieve good numerical convergence at a grid resolution much 
coarser than that required in a fixed, Eulerian grid. In the comoving grid, 
the shock remains at the same location, so the compression rate is applied 
consistently to the CR distribution at the subshock, resulting in much more 
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accurate low energy CR acceleration. 



We showed that the time dependent evolution of CR modified shocks becomes 
approximately self-similar, because the postshock variables including the CR 
pressure approach time-asymptotic values while the shock structure broadens 
linearly with time as /shock ~ Mst/8 and the maximum CR momentum increases 
approximately as pmax ~ khockUs/ i^n- This behavior develops because the ther- 
mal leakage injection rate is controlled mainly by the shock Mach number and 
because the ratio tacc/td is a function of shock Mach number. The self-similar 
behavior is independent of the assumed diffusion coefficient. When the drift 
and dissipation of Alfven waves in the upstream region are included, the CR 
acceleration is less efficient for a given Mach number, as expected from pre- 
vious studies [22|8] . However, at sufficiently high Mach numbers it appears, 
even with the inclusion of the Alfven wave terms, that the efficiency probably 
asymptotes to approximately 50% as measured by both the fraction of energy 
flux transferred to CRs and the fraction of momentum flux that appears as 
postshock CR pressure. 

We summarize here the main properties of our improved simulations including 
Alfven wave effects: 

1. The postshock pressures, Pc and Pg, approach time-asymptotic values quickly. 
At sufficiently strong shocks the postshock CR pressure is about 50 % of the 
shock ram pressure, while the normalized gas pressure drops from its gas- 
dynamic value of 75 %. The normalized postshock gas pressure scales ap- 
proximately inversely with the compression through the shock. The resulting 
total pressure in a strong CR modified shock will be reduced from that in 
an ordinary gas shock of the same initial Mach number. Typically a fraction 
of = 10~^ — 10^'^ of the incoming thermal particles become CRs through 
thermal leakage. 

2. A significant shock precursor develops in response to nonlinear feedback 
from the CR pressure, so that the subshock weakens to Mgubshock ~ 3 for the 
initial shock Mach number Mq > 5. The subshock strength increases very 
slowly with total Mach number in a balance between energy extracted by CRs 
and heating of the precursor. 

3. In our present simulations heating in the precursor is dominated in strong 
shocks by Alfven wave dissipation. The resulting asymptotic density compres- 
sion factor through the precursor can approximated by o"p ~ 1.5(Mo/10)'''^^ 
for Mq > 3, while the compression over the total transition can be approxi- 
mated by at ~ 4.9(Mo/10)'^'^'^ for Mq > 5. This behavior is consistent with 
simple theoretical expectations. Of course, both a.^ and at take gasdynamic 
values for smaller Mq. 

4. Both the CR momentum distribution at the shock, g{xs,p) = fixs,p)p'^, 
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and the integrated distribution. G{p) = J dxg{x,p), clearly exhibit charac- 
teristic concave curvature, reflecting the nonlinear velocity structure in the 
precursor. At non-relativistic momenta G{p) oc with 4.5 < a < 4.2, re- 
flecting the strength of weaker subshocks, while it flattens toward the highest 
momenta. Despite strong nonlinearity in the shock structures, the maximum 
CR momentum still scales reasonably close to predictions from test-particle 
theory. This behavior allows the shock precursor to evolve in an approximately 
self-similar manner. 

5. The CR energy ratio $(Ms) increases with Mq for modest Mach numbers 
and depends on e^, and the ratio 9 = EB/Eth- At sufficiently high Mach 
number $ asymptotes towards ~ 0.5. Since acceleration efficiency is reduced 
by Alfven wave effects, higher Mach numbers are required to approach this 
limit as 9 increases. For larger in the thermal injection model, the mag- 
netic turbulence is weaker, so that particles can stream upstream more easily. 
CR injection is then more efficient. The importance of this effect is significant 
only in relatively weak shocks. The presence of pre-existing CRs is equivalent 
to a higher injection rate. In weak shocks this translates into an apparently 
enhanced acceleration efficiency as measured by the fraction of energy incor- 
porated into CRs. This influence is small in strong shocks, however. 
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Fig. 3. Evolution of a plane-parallel shock with Mq = 10 and T = lO^K 
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number of particles passed through the shock, noUgfit. 
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Fig. 7. Time asymptotic values of the postshock CR and gas pressure in units of the 
shock ram pressure, precursor compression factor, Cp, and total compression factor, 
at, are plotted as functions of the shock Mach number. All models use 9 = 0.1. The 
dotted lines represent Pg 2 and at calculated from the Rankine-Hugoniot relation 
with 7 = 5/3. Symbols are used as follows: open circles for models with Tq = lO^K 
and eB = 0.2, stars for models with Tq = lO^K and = 0.25, and filled triangles 
for models with Tq = lO^k and es = 0.2. 
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Fig. 8. The CR energy ratio ^> is plotted as a function of the shock Mach number 
for different models. In the top panel, es = 0.2 for all models, open circles are 
used for models with Tq = lO^K, while filled triangles are used for models with 
T = lO^K. Long dashed lines are for 9 = 0., solid lines for 6 = 0.1, dotted lines for 
9 = 0.5, and short dashed lines for 9 = 1. cases. In the bottom panel, models with 
open circles and filled triangles are the same as in the top panel and shown only 
to be compared with other cases. Filled circles connected with the solid line are for 
models with Tq = lO^K, €b = 0.2, and 6* = 0.1 with a pre-existing CR population 
with Pcfl/Pgfl = 0.25 — 0.5 and f{p) oc p^^. Stars connected with the solid line: 
To = lO^K, e's = 0.25, and 9 = 0.1. 
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